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A PROGRAM FOR EQUILIBRIUM NORMAL SHOCK AND 

STAGNATION POINT SOLUTIONS FOR 

ARBITRARY GAS MIXTURES 

By Linwood B, Callis and Jane T, Kemper 
Langley Research Center 

SUMMARY 


A computer program written in FORTRAN IV language is presented which yields 
solutions for flow parameters in arbitrary equilibrium gas mixtures in the fol- 
lowing situations: 

(1) Behind a normal shock 

(2) Behind a reflected normal shock 

(3) For in-flight stagnation conditions 

( 4 ) For shock-tube stagnation conditions 

Program output parameters are pressure^ density^ enthalpy^ entropy, compressi- 
bility, temperature, and mole fractions of the included chemical species, which 
may number up to thirty. For traveling normal and reflected shocks, the flow 
velocity and reflected shock velocity are also presented as output. 

Equilibrium flow calculations are carried out by utilizing a free-energy 
minimization technique coupled with the steady-flow conservation equations and 
a modified Newton -Raphson iterative scheme. Chemistry up to second ionization 
is included - 

Input required for the program is described and required physical constants 
for computations involving 27 species of the argon, nitrogen, oxygen, and car- 
bon genre are tabulated. Cases may be run in sequence and any or all of the 
aforementioned flow configurations may be included in a single case. 

Typical shock solutions are presented for argon free air and a model of the 
Mars atmosphere. Air solutions of the present work are compared with those 
from the Avco Corporation and Space Technology Laboratories. 


INTRODUCTION 


Since the advent of hypervelocity vehicles and test facilities^ it has been 
necessary in the solution for normal shock and stagnation point conditions to 
include the effects of the complex chemistry associated with high-speed phenom- 
ena. Since the equilibrium properties of high -temperature air are well known 
(refs. 1 to 6)^ there have been n^Imerous solutions for normal shock parameters 
in air including stagnation point solutions (refs. 7 “to 12 ). Methods of solu- 
tion range from computer use of polynomial fits of equilibrium thermodynamic 
data to hand calculations (for example^ refs. 9 respectively). 

Recently^ however, with the advancing sophistication and success of plane- 
tary probes it has become clear that programs should be available which are 
capable of determining thermochemical equilibrium shock parameters in a gas of 
arbitrary composition. One such program has been developed at the Avco 
Corporation and utilizes a free-energy minimization technique to get the equi- 
librium composition of the flow. This program, however, solves the shock prob- 
lem inversely with values of temperature behind the shock being required in 
order to solve for the remaining shock properties, including shock speed. 

It was believed by the present authors that another useful and somewhat 
more versatile means of solution mi^t be realized by combining free-energy 
minimization techniques with a direct solution to the normal shock problem. The 
RAND method (ref. 13) of equilibrium gas analysis states the equilibrium problem 
simply and, when coupled with steepest descent techniques, is well suited to 
digital-computer use. This method is used in the present program with the con- 
servation equations and a modified Newton-Raphson iterative scheme allowing 
direct solutions for thermodynamic and flow properties behind traveling and 
reflected normal shock waves and at in-fli^t and shock-tube stagnation points. 

The program in FORTRAN IV language and the required program input are 
listed in full detail in appendixes A and B. Typical air solutions are pre- 
sented and comparisons are made with solutions by Ziemer (ref. ll) and Laird 
and Heron (ref. 12 ). In addition, results from normal shock solutions in a 
Martian atmosphere (NASA model 2 , ref. l4) are also presented and compared with 
air solutions. 

The present program is used in conjunction with an IBM 70^*709^ direct 
coupled system. It is referred to as problem 886.5 and is available from the 
Analysis and Computation Division at Langley Research Center, Langley Station, 
Hampton, Virginia. 


SYMBOLS 


fiji,fp perturbation parameters used in stagnation point solution 

H enthalpy 


2 



[(H - T)p] 

(r) 

(n) 

parameter defined by equation (l 6 ) 

L 


characteristic distance 

M 


molecular weight 

Ml 


molecular weight of free-stream gas at 500^ K 

We 


number density of electrons 

P 


pressure 

Ut 


universal gas constant 

R 


xiniversal gas constant divided by free-stream molecular weight 

s/r 


nondimen si onal entropy 

q(k) 

®(n) 


parameter defined by equation ( 5 ) 

[(s - p)h 

(r) 

parameter defined by equation ( 25 ) 

[(s - t)h 

(r) 

parameter defined by equation ( 2 ^]-) 

T 


temperature 

U 


characteristic velocity 

Us^Up 


incident and reflected shock velocity, respectively 

U2 


flow velocity behind traveling incident shock in laboratory 
coordinates 

U2 


velocity relative to shock 

X 


mole fraction 

Z 


compressibility factor, M1/M 

^st 


convergence criterion in stagnation point solution 

P 


density 

PO 


density of standard atmosphere 
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Subscripts: 

1 

2 

5 

5 

st 

sf 

(n) 

P 

i 

J 

Superscripts: 

(k) 

(r) 


character! stic time ^ L/U 
spectroscopic constants 

conditions prior to shock 

conditions behind incident shock in shock tube 
conditions behind standing shock in shock tube 
conditions behind reflected shock 
shock-tube stagnation condition 
in-flight stagnation condition 

refers to minor iteration in shock and stagnation routines 

particular species (atom^ molecule^ ion, and electron) in 
free stream 

particular species (atom, molecule, ion, and electron) 
behind shock 

particular elemental component (atomic elements and 
electrons) 

refers to major iteration in shock routine 
refers to major iteration in stagnation routine 


EQUILIBEIUM PROPERTIES PROGRAM 


The normal shock problem in thermochemical equilibrium is simply solved 
provided there is a straightforward means of handling the necessary equilibrium 
calculations. The RAJNTD method is based on the principle that at given values 
of p and T, a constant-mass equilibriim mixture is so composed that its 
Gibbs free energy is at a minimum value. The composition yielding this minimum 
value is determined by making successive quadratic approximations to the free 
energy and using steepest descent techniques to converge upon a set of mole 
numbers yielding the minimim total free energy. This method permits the formu- 
lation of the equilibriiam problem for an arbitrary gas in direct fashion and 
requires little or no chemical intuition. 
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The equilibrium program used in the present work, requiring p and T as 
thermodynamic input, is a version of the RATOD method and is embodied in work 
done by Allison (ref. 15). This procedure utilizes the free energy minimiza- 
tion technique in conjunction with the partition function of quant-um statistical 
mechanics to determine the free energies and enthalpies of the individual 
species and the equilibrium set of mole numbers (con 5 )osition) . This being 
done, thermodynamic parameters of interest are then determined. Assumptions 
involved in the use of partition functions in this analysis are as follows: 

(1) For molecules the rigid rotor harmonic oscillator model is used, 
account being taken of the variation of vibrational and rotational constants 
due to different electronic configurations. 

( 2 ) Only electronic levels of energy (in the first five electron shells) 
lower than the ionization limit are considered for atoms and atomic ions. 

Effects on thermodynamic properties of vibrational and rotational corrections 
to the model proposed are in general small (approximately 1 percent or less as 
shown in ref. 15 ) and for convenience are neglected in the present work. No 
further details, other than the equilibrium subroutine itself, are presented 
herein on this method of equilibirum gas analysis. Readers interested in these 
details should refer to references I 5 and I 5 . 


ITERATIVE SOLUTION FOR CONSERVATION EQUATIONS 


With an effective equilibrium program available, consideration must now 
be given as to which method of solution of the conservation equations is pref- 
erable, the direct or inverse method. The direct solution requires as input 
the preshock flow conditions, including the shock speed, and yields conditions 
behind the shock. The inverse solution requires the specification of T 
behind the shock with the subsequent solution for the remainder of the condi- 
tions both before and after the shock, including shock speed. The authors of 
the present report believe that the direct method has more general utility. 
Hence, the normal shock problem is approached in this fashion with the aid of 
a modified Newton-Raphson iterative technique (ref. I 6 ) in conjunction with 
the shock-fixed conservation equations. 


Incident Shock 

The conservation equations for the incident shock may be written as 


^ 2^2 = Pi^s 
?2 ^2^ ~ ^1 Pi^s 


( 1 ) 

(2) 
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( 5 ) 


H2 + I u| = Hi + i Uf 

where U 2 is the velocity behind the shock relative to it and values of 
and Ug are specified in the free stream. 

For values of Tj in excess of 800^ values of and Hi, to he used 

in equations (l) to (3)^ are determined with the aid of the equilibrium program. 
Species included in this calculation are those which are to be considered behind 
the shock. This technique makes possible the generation of solutions with dis- 
sociating, high-temperature free streams. For lower values of Tq, the equilib- 
rium program is again used to determine p^_ and Hq; however, the only species 
considered are those present at a temperature of 300^ K. 

Briefly, the iterative solution proceeds as follows: 


(k) (k) 

Pg 

tlon of which is used once again in equations (l) to (3). This pro 


(k) 

(l) Assumed p^ leads, with the aid of equations (l) to (3)^ to 

In principle, these two thermodynamic properties allow the evalua- 

(k+l) 

2 

cedure, which shall be referred to as the major iteration, is repeated to con- 
vergence, solving the problem. In practice, however, the equilibrium program 
requires a p,T input making necessary the Newton -Raphson iterative procedure 
(hereinafter called the minor iteration) in order to determine a value of 
f k^ 

T) ( compatible with the pressure and enthalpy solved for in the major itera- 

vn; 


(k+l) 


tion. With this temperature, or an approximation to it, p 
and the procedure repeated until the desired convergence between 


is determined 
and 


,(k+l) 


is achieved. The superscript (k) refers to the kth major iteration and 


the subscript (n) is associated with the minor iterations. 

( 2 ) The recursive equations used in the minor iteration are 


where 


„(k) „(k) 

= T(k) + ^ ~ M 

(n+l) (n) „(k) 

®(n) 



H 


(k) 

(n) 


T 


(k) 

(n) 


- H 


(k) 

(n-1) 


- T 


(k) 

(n-1) 


( 4 ) 


( 5 ) 
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To begin this iteration, it is necessary to input an estimated teii 5 )erature 
and an initial temperature increment /isT. The sum of these is given as 


T 


( 1 ) 

( 0 ) 


m(l) _ m(l) . Am 
^( 1 ) “ ^( 0 ) ^ 

which is used in the first approximation* Values of the enthalpy H 


( 6 ) 


(k) 

(n) 


in 


equations (h) and (5) are deteniiined from the equilibrium program with 

and as input and are computed with each successive value of as 

(n) (n) 


are 


values of the parameter given by equation (5)- 

(n) 


(3) The program is arranged so that for k = 1 there are three minor 

iterations beginning with the computations of an improved value of 

each successive value of k only one minor iteration is made, use being made 
of the following relations: 


(k+1) ^ g(k) 

(l) last 


(k ^ 1) (7) 


(ktl) ^ (k) 

*^( 1 ) ^last 


(k^l) (8) 


( k) ( k) 

where "^last simply the last values of these quantities com- 

(k) 

puted during the kth major iteration. When the final value of been 


computed. 


m(k) 

(n) 


and 



are used in the equilibrium program to evaluate 


which is compared with to determine whether the specified con- 

vergence criterion has been satisfied. If the criterion has not been satisfied, 
is used in equations (l) to ( 3 ) to obtain and 

cycle is continued until convergence requirements are satisfied. 


Solutions to the problem are presented with reference to a laboratory 
coordinate system. 
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One comment by way of explanation should be made about the scheme of 
minor iterations. The minor iterations are not continued to the point where a 


value of T 


(k) 

(n) 


completely compatible with p^ 


(k) 


and 


H(k) 

2 


is realized. These 


iterations are truncated early in favor of obtaining swiftly a value of p 


(k+l) 


Exp>erience has justified this procedure since it has been obsei*ved that the 
major correction to the estimated temperatuire occurs during the first three 
minor iterations. Subsequent iterations are refinements^ eliminated in favor 
of determining a new approximation to the density and continuing the major 
iteration. This procedure has no bearing on the final accuracy achieved^ inas- 
much as this accuracy is dependent on the convergence criterion specified for 

successive values of . 


Reflected Shock 

The conservation equations for the reflected shock may be written in 
shock-fixed coordinates as 


Up) ~ P^Up (9) 

^ P2^^2 ^ Up) = p^ + P^Up ( 10 ) 

H2 + |(U2 + Urf = H5 + I (ll) 

where all velocities are relative to a laboratory reference system and condi- 
tions behind the incident shock will have been determined from the incident 
shock program. 

For the sake of brevity^ it should suffice to say that the major and minor 
iterative schemes^ including the initial inputs are identical in every detail 
to those for the incident shock. Only the form of the conservation equations 
to be used in the major iteration differs in the two situations. For this 
reason^ no further comment on the reflected shock solution is made except to 
indicate that results are presented in a laboratory reference frame. Such ref- 
erence results in zero flow velocity behind the reflected wave^ as expected. 


In-Flight Stagnation Conditions 

In-flight stagnation conditions are determined by solving the incident 
shock problem and locating a thermodynamic state point so that the conditions 
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( 12 ) 


Hsf = % + I Ug2 


(S/R)sf = (S/r )2 


(15) 


are satisfied. With the present equilibrium program, equations (l2) and (l3) 
are achieved with the greatest facility by means of a two-dimensional Newton- 
Raphson iterative process. This procedure, for purposes of brevity and clarity, 
is presented as a series of comments or computational procedures, as follows: 


(l) The series of iterations is 
and temperatures: 


initiated with the following pressures 


P 


(1) 


= P2 




m(l) 

( 0 ) 



t(i) = 
( 1 ) 



iih) 


where fp and f,p are small perturbation parameters required as input. Super- 
scripts and subscripts within parentheses refer respectively to the number of 
approximations to p^p and the number of iterations on temperature at constant 



( 2 ) By utilizing equations (l4) when applicable and 


where 



((n) ^ 1) ( 15 ) 



( 16 ) 
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I 



((r) > 1) 


(17) 



t(2) 

( 1 ) 


_ Qi(l) 

last 


( 18 ) 


successive values of ten^erature are determined until a value of 
reached so that 


t(i’) 

last 


is 


H 


(r) 

last 


H 


sf 


Hsf 


= ^st 


(19) 


where is a predetermined small quantity. Entropy convergence according 

to 

sf 

is also checked at this point in the solution. If condition (20) is met^ the 
stagnation-point conditions are determined. If convergence is not achieved^ 
use the following proced'ure. 



^ ^st 


( 20 ) 


( 3 ) Determine new values of p 


and T from 


where 



r-’"”] ’ pW.p(-) 


(r ^ 2) 


( 23 ) 
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(S/R) 


(r) 

last 


(S/R) 


(r-1) 

last 


m(^) ^ 

last last 


(r ^ 2) 


(2if) 


(4) With the new values of pressure and temperature reenter step (2) and 
continue the cycle until both enthalpy and entropy meet the convergence test 
(eq,s, (19) and (20))- This procedure, then, satisfies requirements (l2) and 
( 13 ) to the desired degree of accuracy and the problem is complete. 

All values of enthalpy and entropy are determined from the equilibrium 
program. 


Shock-Tube Stagnation Conditions 

The solution for shock-tube stagnation conditions requires tvo incident 
shock solutions followed by a Newton-Raphson iteration in two dimensions. The 
incident shock solution, described previously, yields solutions for conditions 
which serve as input required to determine conditions behind the standing shock. 
This solution, used once again, determines conditions behind the standing shock 
which, when coupled with the two-dimensional iteration described in the previous 
section, are then used to determine a thermodynamic state point so that 


Hst = H2 + 1 u| 

(25) 

(S/R) St = (s/R )5 

(26) 


with the required degree of convergence specified. The solution reached in 
this manner is that of the shock- tube stagnation point. 


The program, briefly described in the preceding sections, is listed in 
FORTRAN IV language in appendix A with a description and explanation of the 
input required given in appendix B. Appendix C presents a compilation of 
physical constants, most of which may be found itemized with regard to their 
source in reference I 5 and references I 7 to 20. These constants are required 
by the program for use with the 27 chemical species indicated. 


LIMITATIONS 


Limitations on the present normal shock program are those restrictions 
placed on the equilibrium properties program. It is recommended in reference I 5 
that the present version of the RAND method of computing the equilibrium com- 
positions be restricted to pressures below 10^ atmospheres; thus, real-gas 
effects manifested at hi^er pressures are avoided. Such effects, for air, are 

11 


I 



taken into account in the work of Lewis and Burgess (ref. 7)^ which allows 
consideration of pressures several orders of magnitude hi^er than the upper 
limit suggested for the present work. 

The present equilibrium program is valid down to pressures at which the 
assumptions Inherent in the theory of statistical thermodynamics begin to fail. 
However^ when a flow process is considered^ pressures must be such that the 
response times of thermodynamic parameters to a change in condition must be 
much smaller than the characteristic time t associated with the problem. 

Such decisions are left to the discretion of the individual Investigator. 

Temperatures considered must be such that only negligible contributions 
are realized from coulomb interactions and from electronic energy levels past 
the fifth electron shelly both these considerations being unaccounted for in 
the equilibrium program. The latter consideration proves to be no problem for 
temperatures below 15 000^ K; however^ this temperature cannot be said to be a 
lower universal limit for the neglect of coulomb effects. These effects 
depend^ in addition to temperature^ on the pressure and the gas mixture con- 
sidered. The problem therefore becomes one for the individual investigator 
depending on the particular circumstances at hand. 


COMPARISON OF SOLUTIONS 


As an indication of the versatility and validity of the present computer 
program^ comparisons are made with the work of Laird and Heron (ref. 12) and 
Ziemer (ref. ll) . 

In figure 1^ the flow configurations used in this study are illustrated. 

In figures 2 to 5 normal shock properties in argon free and carbon free air as 
determined from the present program are compared with those of Laird and Heron 
(ref. 12). The initial mixture considered^ at pj_ = and 10“5 cm Hg and 

T^ = 300^ was composed of 21.153 percent O2 and 78-8^7 percent N2 by volume. 
Species considered were NO*^^ NO^ ^2^ y ^2^ ^2~^ ^ ^2> N“^, N^ 

e”^ N2O excluding N” which was included in reference 12. Figures 2 to 5 indi- 
cate excellent agreement^ small discrepancies being noted in T and p for 
T > 18 000^ K. These discrepancies^ however^ remain below 3 percent. Althou^ 
only data for the incident shock and the shock-tube stagnation conditions are 
shown ^ standing shock and reflected shock data compare equally well. 

In figures 6 and 7^ comparisons are shown between the data of Ziemer 
(ref. 11) and data from the present program in which the previously described 
air model was used. The data of Ziemer was generated by using a graphical rep- 
resentation of the equilibrium air analysis of reference 1. The comparisons^ 
which are made for incident and reflected shock data^ at two initial pressures 
and Tq = 273-2^ K, show considerably larger discrepancies in T and p than 
the Avco data, sometimes reaching 12 percent. These discrepancies, previously 
pointed out by Hoshizaki (ref. 2l),are attributed to errors in Ziemer' s data 
which, according to Ziemer, have an estimated accuracy of 1 to 10 percent. 
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No comparisons^ for identical flow conditions, have "been made of solutions 
discussed in the present section. It was believed that the interpolation 
required between references 11 and 12 would only compound existing 
discrepancies . 

In figure 8, a comparison is made between the incident shock properties of 
the present air model and a Martian atmosphere (NA8A model 2, ref. l4). The 
Mars atmosphere consists of an initial mixture of 10.8 percent CO2 and 89.2 per- 
cent by volume. Species considered in the Martian atmosphere were N, N“*“, 

N++, 0, 0+, 0“*^, 0", C, C+, N 2 , N2'*', O 2 , 02“*", NO, N0+, CO, C0+, CN, 

C02^ e~. Initial conditions of both mixtures were taken to be 
Pj^ = 10“^ atmosphere and = 275.2^ K. 

Finally, comparisons made in reference I 5 of equilibrium air properties 
as generated by the present equilibrixam program with the more rigorous data of 
Gilmore (ref. l) and Browne (ref. 22) generally agree within 1 percent for 
1000° K ^ T ^ 25 000° K and 10"^ ^ p/PO = 10. 


CONCLUDING REIMARKS 


The present program has proved capable of accurately determining both flow 
and thermodynamic parameters behind incident and reflected shock waves and at 
stagnation points for both in-flight and shock-tube thennochemical equilibriixm 
flow. Capable of handling arbitrary free-stream mixtures and gas chemistry 
(up to second ionization), the solution provides a convenient means of solving 
directly for the flow parameters if the free-stream pressure, temperature, 
velocity, and composition are given. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., October 6, I 965 . 
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APPENDIX A 


PROGRAM FOR CALCULATION OF NORMAL SHOCK AND STAGNATION-POINT CONDITIONS 


The program for calculating flow parameters in thermochemlcal equllibrixun 
for nomnal shock and stagnation -point conditions in arbitrary gas mlxtiores was 
written in FORTRAN TV language for the IBM 7^9^ electronic data processing 
system. This program including subroutines and comments is reproduced in the 
following pages. 
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APPENDIX A 


SIBFTC P8865 
C 

C P-886.5 

C NORMAL SHOCK PROGRAM 

C PROGRAMMED FOR THE IBM 7094 

C YIELDING SOLUTIONS FOR FLOW PARAMETERS IN ARBITRARY GAS 

c mixtures in the following SITUATIONS- 

C 1. BEHIND NORMAL SHOCK 

C 2. BEHIND A REFLECTED NORMAL SHOCK 

C 3. FOR IN FLIGHT STAGNATION CONDITIONS 

C 4. FOR SHOCK TUBE STAGNATIONS CONDITIONS 

C 

DIMENSION OMEG (5.30 . 30 ) iF (30 ) *CAPM (30 ) . A ( 1 0 . 30 ) . L ( 30 ) . 

1 G ( 30 » 30 ) »SM ALE (30.30 ) . X ( 30 > » CAPLAM (30.30). DHFO (30 ) . YSTO (30 ) 

C 

D I MENS ION CP (5 ) .BETA (5 ) . AM (5 ) .CODE ( 30 ) . SHBL (8 ) . REBL (8 ) .STBL (8 ) . 
1 STSBL (8 ) .NBTA (5 ) 

DIMENSION TEMP (5) 

DIMENSION YSAVE(30) 

EQUILIBRIUM INPUT 

COMMON OMEG.F.CAPM. A 
COMMON G.SMALE 

COMMON CONH .CONK , CONR . CONPRF . CONNO . N . M . L . N I T 
COMMON EPSl .EPS2. CAPLAM 
COMMON DHFO. AMC. ICl .YSTO 

SHOCK PROGRAM INPUT 

COMMON P10.T10.US.RH02.EPS5.NB.CP.BETA.AM.TP.DELT. IT. NBTA 
REFL INPUT 
COMMON TR.RHOR 


stag input 

COMMON FP.FT.es 

ST-STAG INPUT 

COMMON TSTAG.RHOS 
C 

COMMON SHBL. STBL . STSBL .REBL . SOR . X 
C 

WRITE (6.200) 

200 FORMAT ( IHl .20H EQUILIBRIUM PR0GRAM//27H JANE KEMPER FOR LIN CALLIS 
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1//30H J. O. RGH-125 PROS. NO. 886.5///) 

THE FIRST N+1 data CARDS CONTAIN INFORMATION (IN COLUMNS 1-6) 
identifying each of the N species. the first card is N (COL. 1-3) 


READ(5.205)NCOOE 

205 FORMAT ( 13 ) 

READ (5.206) (CODE( I ) . 1=1 .NCODE) 

206 format ( 1 A6 ) 

C 

C LOAD INPUT All data (AFTER FIRST N+1 CARDS) IS ON DEC CARDS 

C LAST CARD IS FOLLOWED 3Y A TRA 2.4 CARD 

C 

1 CALL LOAD(NERR) 

C 

C NERR CONTAINS LOADING CODE 

C =1 PROPER LOADING. COMPUTE 

C =2 ERROR. CARD IN ERROR WILL BE PRINTED. EXIT 

C =3 END OF FILE. EXIT 

C 

IF (NERR-2 )4 .3.2 

2 CONTINUE 

3 CALL EXIT 

4 N = N 

read (5.205)ITEST 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ITEST CONTAINS CODE FOR OPTIONS (PUNCHED IN COL. 3 ON CARD 
FOLLOWING TRA CARD AFTER LAST DATA CARD) 

POSSIBLE OPTIONS 

0 INCIDENT SHOCK ONLY 

1 INC. AND REFLECTED SHOCKS 

2 INC. AND FREE STREAM STAGNATION 

3 INC.. FREESTREAM STAG.. AND REFL. 

4 INC. AND SHOCK TUBE STAGNATION 

5 INC.. SHOCK TUBE STAG. AnD REFL. 

6 INC.. SHOCK TUBE AND FREESTREAM STAGNATION 

7 INC.. SHOCK TUBE STAG.. FREESTREAM STAG. AND REFL. 


M = M 

CALL SLITE(O) 

RH010=P10/(Tip/300. )*.040619*AMC 
DO 5 I = 1 .NB 
I 1 =N3TA ( I ) 

5 TEMPI I )=CODE( I 1 ) 

WRITE! 6.4000 IRHOIO.PIO.US.TIO. (TEMPI I ). 1 = 1 .NB) 

4000 F0RMAT(1H125H INPUT FOR INCIDENT SH0CK//8H RHO-1 =E15.8.2X. 

15HP-1 =E15.8.2X.3HUS=E15.8.2X.5HT-1 =E15.8//5H BET A// ( 5 ( 9X . 1 A6 . 2X ) 


2 ) ) 

WR I TE (6 .400 1 ) (BETA ( I ) . I = 1 . NB ) 
4001 FORMAT (5E17. 8) 

WR1TE(6.4002)RH02 
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4002 F0RMAT(//13H assumed RhO=E15.§) 

C CONVERSION 

P10=P10*1 .01325E6 
RHOl 0=RHOl 0*1 .E-3 
RH02=RH02*1 .E-3 
US=US*30.48 
DO 3004 I =1 »N 

3004 YSAVE( I )=YSTO( I > 

IF TlO LESS THAN 800 DEGREES KELVIN. COMPUTE ENTHALPY (HIO) 

FROM FREESTREAM COMPOSITION (USING SUBROUTINE ECOM ) 

IF TlO GREATER THAN 800. DEGREES. SOLVE FOR ENTHALPY 
ITERATIVELY USING SUBROUTINE ECOM 

IFCTIO -800.13003.3003.3006 
3003 DO 3005 I =1 .N 

3005 YSTO ( I ) =0. 

SENSE LIGHT 4 USED TO SIGNAL SUBROUTINE ECOM 
C 

CALL SLITE(4) 

DO 3035 1=1 .NB 

J=N3TA ( I ) 

3035 YSTO(J)=BETA( I )/AMC 

3006 C=2.99793E10 

CALL ECOM ( TIO.PIO.OOZ.HOZRT, HIO. RHOl 0 ) 

DO 3007 I = 1 .N 

3007 YSTO ( I) =YSAVE ( I ) 

C 

C STORE INITIAL P.T.RHO.AND U IN SHBL ( 1 -4 ) 

C 

SHBL (1 ) =P1 0 
SHBL (2 ) =T1 0 
SHBL ( 3 ) =RH01 0 
SHBL (4 ) =US 

CALL SHOCK ( TP.OOZ .HOZRT . H2 .SHBL . H 1 0 . NN ) 

C 

C UPON RETURN SHBL ( 5-8 ) CONTAINS P2 . T2 . RH02 . AND UF 

C 

PPRINT=SHBL (5 )/CONPRF 
C 

UPR 1NT= SHBL ( 81/30.48 
C 

WR I TE ( 6 . 4003 lUPRI NT . NN 

4003 FORMAT (////7H 0UTPUT//5H U2 =E 1 5 . 8 . 20X . 25HN0 . OF MAJOR ITERATIONS 
1 = I 3 1 

C 

2050 WR I TE ( 6.202 IPPR I NT . SHBL ( 7 1 . OOZ . HOZRT . SOR . SHBL (6 1 . AMC 
202 FORMAT (//9X. 1 HP . 1 3X . 3HRH0 . 1 4X . 3H 1 /Z . 14X.5HH/ZRT. 12X.3HS/R. 

1 1 4X. 1 HT . 1 6X.2HM1 // (7E1 7.8 1 1 
WRITE(6.221 1 
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221 F0RMAT(//23H FINAL X FROM I TERAT I ON * 4X . 7HSPEC 1 ES// ) 

WRITE (6*220 ) (X ( I ) *CODE ( I ) , I = 1 *N ) 

220 FORMAT(E17.8*10X. 1 A6> 

STORE S/R FOR STAG AND STSTAG SUBROUTINES 

STBL (3 ) =S0R 

TEST FOR ADDITIONAL OPTIONS 

9 IF ( ITEST ) 1 ♦ 1 « 1 0 

10 1C=1 
C 

105 IF ( ITEST-I C ) 1 * 1 1 ♦ 1 2 

11 GO TO ( 15.20* 15*25*25*25*25) ♦ IC 
C 

12 IC=IC+1 
GO TO 105 

C 

15 WRITE(6*4004 )RHOR 

4004 FORMAT( IHl * 16H REFLECTED SH0CK//17H ASSUMED RHO =E15.8) 

RHOR=RHOR*l «E-3 
C 

C STORE P2.T2.RH02, AND UF IN REBL ( 1 -4 ) 

C 

REBL ( 1 ) =SHBL (5 ) 

REBL (2 ) =H2 

REBL (3 ) =SHBL (7 )*1 .E-3 
REBL (4 ) =SHBl (S ) 

CALL REFL ( TR *OOZ* HOZRT *NN ) 

C 

C UPON RETURN REBL(5-8) CONTAINS PR * TR * RHOR * AND UR 

C 

UPRINT = REBL ( 8 )/30 *48 
PPRIN=REBL (5)/C0NPRF 
WR I TE ( 6 * 500 3 ) UPRI NT *NN 

5003 F0RMAT(///7H 0UTPUT//5H U-R= E 1 5 . 8 , 2 OX * 25HN0 . OF MAJOR ITERATIONS = 
113) 

C 

WR I TE (6*202 )PPR1 N *REBL ( 7 ) , 00 Z* HOZRT * SOR* REBL (6 ) , AMC 
WRITE (6*221 ) 

WRITE (6*220 ) (X ( I ) .CODE ( I ) , I = 1 *N) 

ITEST= I TEST-1 
GO TO 10 

20 HS=H10+(SHBL(4)**2)/2. 

WRITE(6*5004 )HS*STBL (3 ) 

5004 FORMAT ( 1 HI . 76H IN FLIGHT STAGNATION POINT DATA IS COMPUTED FROM T 
IHE INCIDENT SHOCK OUTPUT//22H STAGNATION ENTHALPY =E15.8*5X* 
220HSTAGNAT I ON ENTROPY =E15.8/////) 

C 

C STORE P2, AND T2 IN STBL ( 1 -2 ) . H2 IN STBL(4) 
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C 

STBL( 1 )=SHBL(5) 

ST3L(2)=SHBL(6) 

STBL (4 ) =H2 

CALL ST AG ( HS . OOZ « HOZRT . RHO » S TBL « NN ) 

UPON RETURN STBL (5-8) CONTAINS P-ST« T-ST « S/R-ST* AND H-ST 
C 

PPRIN=STBL (5)/C0NPRF 

WR I TE ( 6 » 5005 ) PPR I N « STBL ( 6 ) ♦ NN « OOZ » HOZRT « RHO » S TBL ( 7 ) « STBL ( 8 ) 

5005 FORMAT (///7H OUTPUT//15H PRESSURE =E15.B«18H TEMPERATURE = 

1EI5.8»22X«25HN0« OF MAJOR ITERATIONS = I 3//6H 1/Z =E15.8«5X. 
X7HH/ZRT =E15. 

28,5X.5HRH0 =E 1 5 .8 » 5X . 5HS/R =E1 5« 8 .5X , 3HH =E15.8) 

WRITE (6*221 ) 

WRITE (6 *220 ) (X ( I ) *CODE ( I ) , I =1 *N) 

ITEST=ITEST-2 
GO TO 10 
C 

C TO COMPUTE THE SHOCK TUBE STAGNATION POINT PROPERTIES* SUBROUTINE 

C SHOCK IS USED FOR THE STANDING SHOCK DATA. THEN SUBROUTINE 

C STAG COMPUTES THE STAGNATION POINT DATA. 

C 

C STORE P2*T2*RH02, AND UF IN STSBL ( 1 -4 ) 

C 

25 STSBL ( 1 ) =SHBL (5 ) 

STSBL(2)=SHBL(6) 

STSBL(3)=SHBL(7)*1 .E-3 
RH02=RHOS*T .E-3 
STSBL(4 )=SHBL(8) 

CALL SHOCK (T ST AG* OOZ* HOZRT *HST* STSBL *H2 *NN ) 

HS=H2+(SHBL(8)**2)/2. 

PPRINT=STSBL(5)/C0NPRF 
UPRINT=STSBL(8 )/30.48 
WRITE (6 *5006) 

5006 FORMAT! IHl *84H SHOCK TUBE STAGNATION POINT DATA IS COMPUTED FROM T 
IHE FOLLOWING STANDING SHOCK DATA) 

WR 1 TE (6*202 ) PPR I NT * STSBL ( 7 ) * OOZ * HOZRT * SOR * STSBL ( 6 ) * AMC 
WRITE(6*5007)HS 

5007 F0RMAT(//39H SHOCK TUBE STAGNATION POINT ENTHALPY =EI7.8/) 

C 

C STORE P*T* AND H VALUES FROM SECOND ENTRY INTO INCIDENT SHOCK 

C SUBROUTINE IN STSBL (1-2*4) 

C 

C STORE S/R-2 IN STSBL (3) 

C 

STSBL! 1 )=STSBL(5) 

STSBL ( 2 ) =STSBL ( 6 ) 

STSBL(3)=SOR 
STSBL(4 )=HST 

call stag ( HS * OOZ * HOZRT * RHO * STSBL * NN ) 
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C 

UPON RETURN STSBL(5-8) CONTAINS P-STS. T-STS* S/R-STS» AND H-STS 
PPRINT=STSBL (5 )/CONPRF 

WRlTE(6t5005)PPRINT.STSBL(6 ) « NN» OOZ ♦ HOZRT . RHO . STSBL ( 7 ) « STSBL ( 8 ) 
WRITE (6 <221 ) 

WRITE (6 <220 ) (X( 1 ) «CODE ( I ) « 1=1 *N) 

ITEST= I TEST-4 
GO TO 10 
END 


IBFTC SHOCK 

SUBROUT I NE SHOCK ( TGUESS « OOZ . HOZRT . H2 * BLOCK , H 1 0 ♦ NCOUNT ) 

THIS SUBROUTINE USES A ONE-DIMENSIONAL NEWTON-RAPHSON ITERATION 
SCHEME TO FIND TEMPERATURE AND PRESSURE AT EQUILIBRIUM BEHIND 
INCIDENT SHOCK. IT WILL CALL SUBROUTINE ECOM TO COMPUTE THE 
EQUILIBRIUM PROPERTIES. 

HIO IS INITIAL ENTHALPY 
TGUESS IS TEMPERATURE ESTIMATE 

BLOCK (1-4) contains INITIAL P.T. DENSITY AND VELOCITY 

FINAL values OF P.T.DENSTIY AND VELOCITY STORED IN BLOCK (5-8) 
ENTHALPY STORED AT H2 
NCOUNT IS an iteration COUNT 
1/Z. H/ZRT STORED IN OOZ AND HOZRT 


DIMENSION OMEG (5. 30 . 30 ) . F ( 30 ) ♦ CAPM ( 30 ) . A ( 1 0 . 30 ) . L ( 30 ) . 

1 G ( 30 . 30 ) . SMALE ( 30 ♦ 30 ) « X ( 30 ) . CAPLAM (30.30). DHFO ( 30 ) . YSTO ( 30 ) 
DIMENSION CP(5 ) .BETA (5 ) . AM (5 ) .SHBL (8) .REBL (8 ) .STBL (8 ) .STSBL (8 ) . 
INBTA (5 ) 

DIMENSION T(2) .H(2) .BLOCK(8) 

COMMON OMEG.F.CAPM. A 
COMMON G. SMALE 

COMMON CONH . CONK . CONR . CONPRF . CONNO . N . M . L . N I T 
COMMON EPSl .EPS2, CAPLAM 
COMMON DHFO . AMC . 1 C 1 . YSTO 

COMMON P 1 0 . T 1 0 . US . RH02 . EPS5 . NB . CP . BETA , AM . TP . DELT . I T . NBT A 
COMMON TR.RHOR 
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C6mmON FPtFT«eS 
COMMON TSTAG*RH0S 
C 

COMMON SHBL * STBL « STSBL » REBL * SOR t X 
LET RH02 = FIRST RHO 
VELl (AA )=C*D/AA 

PRESl (AA,BB )=B+C*D**2-AA*BB**2 
ENTHl ( AA ) =H1 0+ (0**2 )/2.- ( AA**2 )/2. 
NCOUNT= 1 
B=BL0CK ( 1 ) 

C=BL0CK<3) 

D=BLOCK (4 ) 

U2 = VEL1 (RH02 > 

P2=PRES1 (RH02*U2) 

H2=ENTH1 (U2 ) 

COMPUTE FIRST POINT 

ITT=IT 

7 T ( 1 ) =TGUESS 

CALL ECOM ( T ( 1 ) » P2 «00Z i HOZRT , H < 1 ) »RHO ) 
COMPUTE SECOND POINT 
T(2)=T( 1 )+DELT 

CALL EC0M(T(2) , P2 . OOZ » HOZRT ♦ H ( 2 ) *RHO) 
S=(H(2)-H(1 ))/(T(2)-T(l )) 

T ( 1 ) =T (2 ) 

TEMPERATURE FROM FIRST ITERATION 

T ( 2 ) =T (2 )+ (H2-H (2 ) ) /S 
H ( 1 ) =H ( 2 ) 

IF(T(2) )25»25»8 

8 CALL EC0M(T(2) »P2»00Z*H0ZRT,H(2) »RHO) 

S IS SLOPE (H2-H1 )/'(T2-Tl ) 

85 S= (H(2)-H( 1 ) )/(T(2)-T( 1 ) ) 

T3 

T( 1 >=T(2) 

TEMPERATURE FROM SECOND ITERATION 

T(2)=T(2)+(H2-H(2) ) /S 
H ( 1 ) =H (2 ) 

IF(T(2 ) )25*25»9 

9 CALL EC0M(T(2> «P2*00Z»H0ZRT.H(2) *RHO) 
C 
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C IF ITT IS GREATER THAN 2* ITERATE AGAIN ON TEMPERATURE WITH 

C FIRST PRESSURE 

C 

IF( ITT-2)10* 10. 1 1 
10 SLAST=(H(2)-H(1 ) )/(T(2)-T(1 ) ) 

TLAST=T (2 ) 

GO TO 12 

c 

1 1 ITT=ITT-1 
GO TO 85 
C 

c 

C TEST RHO FOR CONVERGENCE 

C 

12 IF(ABS( (RH0-RH02)/RH02 )~EPS5 >20^20# 13 
C 

C NON-CON VERGENCE- 

C 

C COMPUTE NEW PRESSURE AND CONTINUE ITERATION ON TEMPERATURE AND 

C PRESSURE UNTIL RHO CONVERGES 

C 

13 RH02=RH0 
U2=VEL1 (RH02 ) 

P2=PRESI (RH02.U2) 

H2=ENTH1 (U2 ) 

NC0UNT=NC0UNT+1 

T { 1 ) =TLAST 

14 CALL ECOM ( T ( 1 ) .P2 »00Z. HOZRT . H( 1 ) .RHO ) 

145 S=SLAST 

15 T <2 ) =T ( 1 )+ (H2-H ( 1 ) )/S 
CALL ECOM (T (2) . P2 . OOZ . HOZRT ♦ H ( 2 ) .RHO) 

SLAST= (H(2)-H(1 ))/(T(2)-T(1 ) ) 

TLAST=T (2 ) 

GO TO 12 

CONVERGENCE - STORE OUTPUT 

20 U2=VELI (RHO ) 

UF=US-U2 
rH0=RH0*1 .E3 
BL0CKC5 )=P2 
BLOCK (6 ) =TLAST 
BLOCK(7)=RHO 
BLOCK (8 ) =UF 
RETURN 

TEMPERATURE ESTIMATE TOO HIGH - ADJUST 

C 

25 TGUESS= (TGUESS-T 1 0 ) /2. 

GO TO 7 
END 
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SIBFTC REFL 

SUBROUTINE REFL ( TGUESS ♦ OOZ ♦ HOZRT * NCOUNT ) 

C 

C THIS SUBROUTINE USES A ONE-DIMENSIONAL NE WTON-RAPHSON ITERATION 

C SCHEME TO FIND TEMPERATURE AND PRESSURE At EQUILIBRIUM BEHIND 

C REFLECTED SHOCK. IT WILL CALL SUBROUTINE ECOM TO COMPUTE THE 

C EQUILIBRIUM PROPERT lES# 

C 

DIMENSION OMEG(5,30.30) ♦ F ( 30 ) i CAPM ( 30 ) » A ( 1 0 « 30 ) . L ( 30 ) i 
1 G (30.3 0 ) .SMALE ( 30 . 30 ) . X ( 30 ) .CAPLAM (30.30 ) .OHFO (30 ) . YSTO (30 ) 

D I MENS ION CP (5 ) .BETA (5 ) . AM ( 5 > . SHBL (8 ) .REBL (8 ) .STBL (8 ) . STS0L (8 ) . 
INBTA (5 ) 

DIMENSION T(2) .H(2) 

C 

COMMON OMEG.F.CAPM. A 
COMMON G.SMALE 

COMMON CONH . CONK , CONR . CONPRF ♦ CONNO . N . M . L . N I T 
COMMON EPSl .EPS2.CAPLAM 
COMMON DHFO.AMC. ICl .YSTO 

COMMON PI 0 . T 10 .US .RH02 .EPS5. NB. CP .BETA. AM. TP. DELT . I T.NBTA 
COMMON TR . RHOR 
COMMON FP.FT.ES 
COMMON TSTAG.RHOS 
C 

COMMON SHBL . STBL . STSBL . REBL . SOR . X 
C 

C LET RHOR = FIRST RHO 

C 

VEL ( A A ) =C*D/ ( AA-C) 

PRES ( AA .BB ) =B+C* ( AA+D ) **2-BB*AA**2 
ENTH ( AA ) =E+.5* ( D+AA ) **2-.5*AA**2 
NCOUNT= 1 
0=REBL ( 1 ) 

C=REBL ( 3 ) 

D=REBL(4 ) 

E=REBL (2 ) 

UR=VEL (RHOR ) 

PR=PRES (UR. RHOR ) 

HR=ENTH (UR ) 

COMPUTE FIRST POINT 
C 
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7 T( 1 )=TGUESS 

CALL ECOM(T( 1 ) , PR « OOZ ♦ HOZRT t H ( 1 ) ♦RHO) 

C 

C COMPUTE SECOND POINT 

C 

T(2)=T(1 )+DELT 

CALL EC0M(T(2) , PR t OOZ ♦ HOZRT , H ( 2 ) »RHO) 

C 

C S IS SLOPE (H2-H1 )/(T2~Tl ) 

C 

S=(H(2)-H(1 ))/{T(2)->T(l)) 

T ( 1 ) =T (2 ) 

C 

C TEMPERATURE FROM FIRST ITERATION 

C 

T (2 ) =T (2 )+ (HR-H (2 ) ) /S 
H{ 1 ) =H (2 ) 

IF(T(2) )25f 25«8 

8 CALL EC0M(T(2) fPR500Z<H0ZRT,H<2 ) fRHO) 

85 S= (H(2)-H( 1 ) )/(T(2)“T( 1 ) ) 

T ( 1 ) =T (2 ) 

C 

C TEMPERATURE FROM SECOND ITERATION 

C 

T (2 ) =T (2 )+ (HR-H (2 ) )/S 
H(1 )=H(2) 

IF(T(2) )25*25»9 

9 CALL ECOM ( T (2 ) ^PR«OOZ*HOZRT f H(2 ) fRHO ) 

C 

C IF ITT IS GREATER THAN 2< ITERATE AGAIN ON TEMPERATURE WITH 

C FIRST PRESSURE 

C 

IF ( IT-2 ) 10 t 1 0 f 1 1 

10 SLAST= ( H (2 ) -H ( 1 ) )/ ( T (2 )-T ( 1 ) ) 

TLAST = T (2.) 

GO TO 12 
C 

1 1 ITT=ITT-1 
GO TO 85 

C 

C TEST RHO FOR CONVERGENCE 

C 

12 IF ( ABS ( (RHO-RHOR ) /RHOR ) -EPS5 )20« 20* 1 3 
C 

C 

C NCN-CONVERGENCE- 

C COMPUTE NEW PRESSURE AND CONTINUE ITERATION ON TEMPERATURE AND 

C PRESSURE UNTIL RHO CONVERGES 

C 

13 RHOR=RHO 
UR=VEL (RHOR ) 
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PRs:PRE:5 (UPi phop i 

HR=ENTH (UR ) 

NCOUNT=NCOUNT+l 
T( 1 )=TLAST 

CALL ECOM(T( 1 ) ♦ PR f OOZ # HOZRT ♦ H ( 1 ) fRHO) 

S=SLAST 

T (2 ) =T (2 ) + ( HR-H (2 ) ) /S 

CALL EC0M(T(2) , PR , OOZ * H02RT • H ( 2 ) «RHO) 

SLAST= (H (2 )-H(l ))/CT(2)-T(l ) ) 

TLAST=T (2 ) 

GO TO 12 
C 
C 

CONVERGENCE - STORE OUTPUT 

20 UR=VEL ( RHO ) 

RHO=RHO*I .E3 
REBL (5 ) =PR 
REBL (6 ) =TLAST 
REBL (7)=RH0 
REBL (8 ) =UR 
RETURN 
C 

C TEMPERATURE ESTIMATE TOO HIGH - ADJUST 

C 

25 TGUESS= ( T GUESS- T1 0)/2. 

GO TO 7 
END 


SIBETC STAG 

SUBROUT I NE STAG CHS *OOZ « HOZRT • RHO # BLOCK ♦ NN ) 

C 

C THIS SUBROUTINE USES A TWO DIMENSION NEWTON-RAPHSON ITERATION 

C SCHEME TO COMPUTE FREESTREAM OR SHOCK TUBE STAGNATION POINT 

C PRESSURE AND TEMPERATURE GIVEN STAGNATION ENTROPY AND ENTHALPY. 

C THIS ROUTINE CALLS SUBROUTINE ECOM TO COMPUTE EQUILIBRIUM PRO- 

C PERTIES. 

C 

C HS CONTAINS STAGNATION ENTHALPY 

C BLOCKS 1-4) CONTAINS INITIAL VALUES OF P»TfS/Rf AND H 

C FINAL P.T.S/R, AND H ARE STORED IN BL0CK(5-8) 

C NN IS an ITERATION COUNT 

C 1/Z» H/2RT AND DENSITY STORED IN OOZ« H0ZRT« AND RHO 
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C 


C 


c 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 


DIMENSION OMEG (5« 30 » 30 ) «F(30 ) »CAPM ( 30 > ♦ A ( 1 0 » 30 ) »L ( 30 ) * 

1G<30«30 ) *SMAl-E ( 30 * 30 ) « X ( 30 ) » CAPLAM (30«30 ) . DHFO (30 ) » YSTO (30 ) 

DIMENSION CP (5 ) «BETA (5 ) « AM (5 ) t CODE (30 ) «SHBL (8 ) tREBL (8 ) «STBL (8 ) ♦ 
1 STSBL (8 ) »NBTA (5 ) 


DIMENSION T(2)«H(2).P(2) «SR(2) tBLOCK (8 ) 

COMMON 0MEG.F«CAPM« A 
COMMON G«SMALE 

COMMON CONH , CONK , CONR « CONPRF t CONNO « N « M . L « N I T 
COMMON EPSl «EPS2,CAPLAM 
COMMON DHFO « AMC « I C 1 . YSTO 

COMMON P I 0 ♦ T 1 0 « US »RH02 « EPS5 ♦ NB « CP « BETA ♦ AM « TP « DEL T 
COMMON TR.RHOR 
COMMON FP,FT»ES 
COMMON TSTAG«RHOS 


TtNBTA 


COMMON SHBL . STBL , STSBL « REBL » SOR . X 


IR = 0 
NN=1 


scheme 

INSIDE LOOP- 
OUTSIDE LOOP- 


ITERATE ON T UNTIL 
ITERATE ON P UNTIL 


ENTHALPY (H) 
entropy (S/R 


CONVERGES 

CONVERGES 


FIRST PRESSURE CURVE 


5 


S0RS=BL0CK(3) 

P(1 )=BLOCK( 1 ) 

T ( 1 ) =BLOCK (2 ) 

SORS=BLOCK (3) 

H( 1 )=BLOCK ( 4 ) 

T (2 ) =T ( 1 )+T ( 1 )*FT 

CALL ECOM (T (2 ) .P( 1 ) ♦ OOZ » HOZRT ♦ H ( 2 ) 
S=(H(2)-H(1 ))/(T(2)-T(l)) 

T ( 1 ) =T (2 ) 

T (2 ) =T ( 1 )+ (HS-H (2 ) )/S 
H(1 )=H(2) 

CALL ECOM (T (2) ,P( 1 ) , OOZ . HOZRT * H ( 2 ) 
S=(H(2)-H(1 ))/(T(2)-T(l)> 

T ( 1 ) =T (2 ) 

H( 1 ) =H (2 ) 


,RHO > 


.RHO) 


REL=ABS (H (2 )-HS )/HS 
lF(REL-ES)15.15ilO 
15 REL=ABS ( (SORS-SOR )/SORS ) 
IF(REL-ES) 1 00. 100 1 20 
10 T(2)=T(1 )+(HS-H(2) )/S 
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66 TO S 

SECOND PRESSURE CURVE 

20 SR ( 1 )=SOR 
SLAST=S 
TLAST = T (2 ) 

P(2)=P( 1 )+P( 1 )*FP 
T(1 )=BLOCK(2) 

23 T (2)=T ( 1 )+ (HS-H ( 1 ) )/SLAST 
TWO DIMENSIONAL ITERATION LOOP 

24 CALL EC0M(T(2) «P(2) ♦00ZfH0ZRTfH(2) ♦RHO) 

IF ( IR)241 ,241 *242 

241 SLAST=(H<2)-H(1 ) )/(T(2)-T(1 ) ) 

242 H( 1) =H (2 ) 

IR=0 

REL = A0S ( (HS-H <2 ) )/HS ) 

IF (REL-ES )30f30 ^25 
C 

C H NON-CONVERGENT» continue ITERATION ON T 

C 

25 T(1 )=T(2) 

GO TO 23 

C 

C H CONVERGES • TEST S/R 

C 

30 SR(2)=SOR 

REL=ABS( (S0RS-SR(2) )/S0RS) 

IF(REL-ES ) lOOf I00f35 
C 

C S/R NON-CONVERGENT - ADJUST P AND T AND CONTINUE ITERATIONS 

C 

35 SP=(SR(2)-SR(1 ) )/(P(2)-P(l ) ) 

ST= (SR (2 )-SR ( I ) )/ (T (2 )-TLAST ) 

P(1 )=P(2) 

TLAST = T (2 ) 

T ( 1 ) =T (2 ) 

P(2)=P( 1)+ (SORS-SR{2 ) )/SP 
T(2)=T(2) + (SORS-SR(2 ) )/ST 
SR(1 )=SR(2) 

NN=NN+1 
IR = 1 

GO TO 24 
C 
C 

C S/R CONVERGES - STORE OUTPUT 

C 

100 BL0CK(5)=P(2) 

BLOCK(6)=T(2) 
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C 

C 


BL0CK(7)=SR(2) 

BL0CK(8)=H<2) 

RHO=RHO*l .E3 

RETURN 

END 


SIBFTC ECOM 

SUBROUT I NE ECOM ( T tPSTO , OOZ « HOZRT ♦ H « RHO ) 

C 

C SUBROUTINE WHICH, GIVEN A TEMPERATURE AND PRESSURE, COMPUTES 

C THE thermodynamic EQUILIBRIUM PROPERTIES OF A GAS DESCRIBED BY 

C THE INPUT. 

C 

C T - TEMPERATURE 

C PSTO - PRESSURE 

C OOZ - 1/Z 

C HOZRT - H/ZRT 

C H - ENTHALPY 

C RHO - DENSITY 

C 

C MOLE FRACTIONS (X(I)) STORED IN COMMON 

C 

DIMENSION OMEG( 5,30,30 ) ,F(30 ) ,CAPM (30 ) , A ( 1 0 » 30 ) , L ( 30 ) , 

IG <30,30 ) ,SMALE (30 ,30 ) , X ( 30 ) , 

2E (30 ) , Y (30 ) ,0 (30) , CAPLAM ( 30 , 30 ) , CAPF I (30>,R(10,10),B<10)» 

3TEMPS (10) ,BSUM (11,1) , ABLOCK (11,11) ,PTEMP (30 ) , ZETA ( 30 ) » 

4ZETAPR ( 30 ) ♦ DHFO ( 30 ) , AL AM ( 30 ) , 

5YST0 (30 ) , IPl VOT (11), INDEX (11 , 2 ) , DQ I NT ( 30 ) , Q I NT ( 30 , 30 ) 

DIMENSION CP (5 ) .BETA (5 ) , AM ( 5 ) ,SHBL ( 3 ) , REBL ( 8 ) , STBL ( 8 ) ,STSBL (8 ) , 
1NBTa'(5 ) 

C 

COMMON OMEG.F.CAPM, A 
COMMON G.SMALE 
COMMON CONH.CONK.CONR 

COMMON CONPRF , CONNO , N , M , L , NIT 

COMMON EPSl ,EPS2,CAPLAM 
COMMON DHFO.AMC, ICl ,YSTO 

COMMON P10,T10,US,RH02,EPS5,NB,CP,BETA, AM.TP.DELT, IT.NBTA 
COMMON TR.RHOR 
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Common 

COMMON TSTAG^RHOS 
C 

COMMON SHBL t STBL » STSBL ♦ REBL « SOR * X 
C 

PI=3. 14159 
C=2.99793E1 0 
NC0UNT=0 
LTEST=LTFST 
N2=N 

DO 5 I = 1 »N 
5 Y( I )=YSTO( I ) 

P=PSTO 

34 TK=CONK*T 
RT=CONR*T 

346 YBAR=0.0 

DO 347 1=1 «N 

347 YBAR=YBAR+Y ( I ) 

DO 40 I =1 ,N 
TEMPI =0 

LEND = L ( I ) 

DO 37 Ll=l fLEND 
IF(F( I ) )31 ,35*31 

31 PROD=l« 

DO 33 IC=1 , ICl 

IF(OMEG( IC,L1 ♦ I ) )32,33,32 

32 PROD=PROD* ( 1 .-EXP (-CONH*C*OMEG ( I C ,L 1 • I ) /TK ) ) 

33 CONTINUE 

PART= ( T/ (CAPLAM (Ll * I > *PROD ) > **F ( I ) 

GO TO 36 

35 PART=1 • 

36 QINTCLl , I )=PART*G(L1 * I ) *EXP ( -CONH*C*SMALE ( L 1 » I )/TK) 

37 TEMP1=TEMP1+QINT(L1 , I ) 

Q ( I ) = ( SORT (2. *P I /CONH*TK/ ( CONH*CONNO ) *CAPM ( I ) )**3 ) *TK/CONPRF* TEMP I 
IF(Y( I )/YBAR >38,38, 39 

38 CAPFI ( I > =0 
GO TO 40 

39 CAPE I ( I ) =Y ( I )* ( AlOG ( P/CONPRF )+ALOG ( Y ( I )/YBAR )-ALOG ( Q ( I ) )+DHFO ( I ) 

1 /RT ) 

40 CONTINUE 
C 

C SENSE LIGHT 4 ON DO NOT ITERATE 

C SENSE LIGHT 4 OFF - ITERATE 

C 

CALL SLITET(4,JJ) 

GO TO (95,396),JJ 
396 DO 50 J=1 ,M 
DO 50 K=1 ,M 
R(K, J)=0.0 
B( J)=0.0 
DO 50 I =1 ,N 
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B(J)=B(J)+A(J* I )*Y( I ) 

50 R(Kf J)=R(K* J)+A ( J* I )*A (Kf I )*Y( I ) 

C 

C SET UP matrix for SOLUTION OF EQUATIONS 

C 

DO 60 J=1 ,M 
TEMPS (J )=0#0 
DO 55 1=1 *N 

55 TEMPS( J)=TEMPS( J)+A( Jf I )*CAPFI ( I ) 

BSUMCJ* 1 )=B( J)+TEMPS(J) 

C 

C constant terms in BSUM BLOCK 

C 

DO 56 K=1 fM 
K1 =K+1 

56 ABLOCK( J*K1 )=R(K» J) 

C 

C PI TERMS IN ABlOCK IN COLUMNS 2 THROUGH N+ 1 

C 

60 ABL0CK( J,1 )=B(J) 

C 

C (X/Y) TERMS IN FIRST COLUMN 

C 

Ml =M+1 

ABL0CK(M1 , 1 )=0.0 
DO 61 K=1 tMl 
K1 =K+1 

61 ABLOCK (Ml ,K1 )=B(K ) 

BSUM (Ml ♦ 1 ) =0*0 

DO 62 I =1 f N 

62 BSUM (Ml ♦ 1 )=BSUM (Ml , 1 )+CAPFl ( I ) 

C 

C MATINV EXPECTS AN M+1 BY M+1 MATRIX 

C 

CALL MAT I NV (ABLOCK ( 1 * 1 ) ^Ml t BSUM ( 1 , 1 ) • 1 tDETERM* IPI VOT< INDEX» 1 1 *0 ) 
C 

C RETURN WITH ANSWERS IN BSUM 

C 

ZETAP=BSUM( 1 f 1 )*YBAR 
ZERO=0% 

NEG=0*0 
DO 70 I =1 »N 
PTEMP( I )=0#0 
DO 65 J=1 f M 
J1 =J+1 

65 PTEMP( I )=PTEMP( I )+BSUM( J1 , 1 ) *A ( J » I ) *Y (I ) 

ZETA ( I ) =-CAPFI ( I )+Y ( I )*35UM (1*1 )+PTEMP ( I ) 

TEST FOR NEGATIVE OR ZERO ZETA 
C 

68 IF (ZETA ( I ) ) 69 *695 *70 
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6 ^ P1ECE=-Y( I )/(ZETA(I)-Y< I ) ) 

IF (PIECE ) 691 • 692 « 691 

691 NEG=NEG+1 
ALAM(NEG)=PIECE 
GO TO 70 

692 Y( I )=0 
ZERO=l # 

GO TO 70 

69*5 IF(Y( I ) >69,70 #69 

70 CONTINUE 
C 

C FIND GREATEST NEGATIVE ZETA-Y 

C 

I F ( ZERO ) 700 ♦ 700 ♦ 698 

698 I F (NCOUNT-N I T ) 699 ,100,100 

699 NCOUNT=NCOUNT+l 
GO TO 346 

700 IF(NEG-1 >78,71 ,73 

71 ALAmPR= #999999*ALAM ( 1 > 

GO TO 745 

73 ARGl =ALAM ( 1 > 

DO 74 I=2,NEG 

72 ARG2=ALAM ( I ) 

ARGl = AM INI ( ARGl , ARG2 > 

74 CONTINUE 

AlAMPR= •999999*ARG1 
745 I IC=0 

75 ZETAP=0 

DO 76 I =1 ,N 

ZETAPR( I >=Y( I )+AlAMPR* (ZETA ( I > ~Y ( I > > 

76 ZETAP = ZETAP + ZETAPR ( I ) 

OLAM=0 

DO 77 I =1 ,N 

IF(ZETAPR( I >/ZETAP>77,77,765 

765 DLAM = DLAM+ ( ZETA ( I > - Y ( I > > * ( ALOG (P/CONPRF ) -ALOG ( Q ( I > > +DHFO ( I > /RT + ALO 
1 G(ZETAPR ( I )/ZETAP ) ) 

77 CONTINUE 

IF (DLAM >81 ,81 ,80 

80 IF ( I I C-3 >805,81 ,81 
805 I IC= I I C+1 

ALAMPR=ALAMPR*.9 
GO TO 75 

78 ALAMPR=1% 

GO TO 745 

CONVERGENCE TEST FOR Y(I)S 

81 IF ( ALAMPR-# 50 >83,815,815 
815 DO 82 I =1 ,N 

IF(ZETAPR( I >>8l3,8l6*ai3 
813 REL=Y ( I >-ZETAPR ( I > 
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IF{ABS(REL)~EPS1 )818«ei8«83 
818 REL=ZETAPR( I )/Y( I )-l . 

I F ( ABS { REL ) -EPS2 ) 82 . 82 » 83 

816 IF(Y( I ) )817.82,817 

817 GO TO 83 

82 CONTINUE 
C 

C Y(I)S CONVERGE 

C 

DO 800 1=1 ,N 

800 Y( I )=2ETAPR( I ) 

GO TO 95 
C 

C NON-CONVERGENCE OF Y(I)S 

C 

83 NCOUNT=NCOUNT+l 
IF(NCOUNT-NIT)84« 100.100 

84 DO 85 r=l,N 

85 Y( I )=ZETAPR( I ) 

C 

C REPEAT WITH NEW Y(I)S AND NO. OF ITERATIONS LESS THAN NIT 

C 

GO TO 346 
95 DO 20 1 1 = 1 , N 

201 X ( I ) =Y ( I )*CAPM ( I ) 

YBAR=0.0 
CAPMI =0 
DO 2026 1=1 .N 
YBAR = YBAR+Y ( I ) 

2026 CAPMI=CAPMI+X( I )/CAPM( I ) 

CAPMI=1 .0/CAPMI 
Z=AMC/CAPMI 

ESUM = 0 

DO 2029 I =1 .N 
OSUM=0 
DQINT( I )=0 
LEND=L ( I ) 

DO 2028 Ll=l .LEND 
SUM = 0 

DO 2027 IC= 1 . 1 Cl 

HOOTK=CONH*C*OMEG { IC .L 1 . I )/TK 
IF (OMEG ( IC.Ll . I ) >2000.2027.2000 
2000 SUM=SUM+HOOTK/ (EXPIHOOTK )-l . ) 

2027 CONTINUE 

DO I NT ( I )=DQINT( I ) +Q I NT (L I . I ) * ( F ( I ) XT* ( 1 .+SUM)+SMALE (LI . I ) *CONH*C 
1 X (TK*T ) ) 

2028 QSUM=QSUM+QINT (LI . I ) 

E ( I ) = 1 ./CAPM ( I )* ( 1 ,5*RT + RT»T/QSUM*DQINT ( I )+DHFO ( I ) ) 

2029 ESUM=ESUM+X ( I )*E ( 1 ) 

HOZRT=CAPMI *ESUM/RT+1 .0 
H=HOZRT*CONR*T*Z/AMC 
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TK = t*dON«; 

FSUM=0 

00 2040 1=1 «N 

2033 IF(Y( I )/YBAR)2034«2034«2035 

2034 CAPFI ( 1 )=0 
GO TO 2040 

2035 CAPFI ( 1 )=Y( I ) * ( ALOG ( P/CONPRF )+ALOG ( Y < I ) /YBAR ) - ALOG < Q ( I ) )+DHFO(I ) 
1/RT) 

2040 FSUM=FSUM+CAPFI ( I ) 

SOZR=HOZRT-CAPM 1 *FSUM 

SOR=SOZR*Z 

RHO=P*CAPMI/RT 

U=CAPX+.43429*AL0G(273. 16/(Z*T) > 

OOZ=l .0/Z 
DO 300 I =1 . N 

300 X( I )=X( I )*CAPMI/CAPM ( I ) 

RETURN 

100 WRITE(6.5000) 

5000 FORMAT ( 1 HO « 25H THIS CASE NON-CONVERGENT ) 

CALL EXIT 
END 


SIBMAP 

* 

* 

* 

SYMBOL 


SYMBOL 150 

SUBROUTINE WHICH DEFINES 
LOADING ROUTINE (PW-LOAD 


ENTRY 

SYMBOL 

BCI 

1 fOMEG 

PZE 

OMEG 

BCI 

1 fF 

PZE 

F 

BCI 

1 ♦ capm 

PZE 

CAPM 

BCI 

1 • A 

PZE 

A 

BCI 

1 »G 

PZE 

G 

BCI 

1 *smale 

PZE 

SMALE 

BCI 

1 ♦CONH 

PZE 

CONH 


NPUT AND ITS STORAGE FOR THE 
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0CI 

1 4 CONK 

PTE 

CONK 

BCI 

1 ♦ CONR 

PTE 

CONR 

BCI 

1 ♦ CONPRF 

PTE 

CONPRF 

BCI 

1 f CONNO 

PTE 

CONNO 

BCI 

1 f N 

PTE 

N 

BCI 

1 f M 

P2E 

M 

BCI 

1 *L 

PZE 

L 

BCI 

1 f NIT 

PZE 

NIT 

BCI 

1 ^EPSl 

PZE 

EPSl 

BCI 

1 f EPS2 

PZE 

EPS 2 

BC I 

1 ♦ CAPLAM 

PZE 

CAPLAM 

BCI 

1 ♦ DHFO 

PZE 

DHFO 

BCI 

1 f AMC 

PZE 

AMC 

BCI 

1 • ICl 

PZE 

ICl 

BCI 

1 » YSTO 

PZE 

YSTO 

BCI 

1 #P10 

PZE 

PI 0 

BCI 

1 f TIO 

PZE 

TIO 

BCI 

C 

(/) 

PZE 

US 

BCI 

1 fRH02 

PZE 

RH02 

BCI 

1 •EPSS 

PZE 

EPS5 

BCI 

1 fNB 

PZE 

NB 

BCI 

1 fCP 

PZE 

CP 

BCI 

1 fBETA 

PZE 

BETA 

BCI 

1 • AM 

PZE 

AM 

BCI 

1 f TP 

PZE 

TP 

BCI 

I f DELT 
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d^Lt 


BCI 

1 • IT 


PZE 

IT 


BCI 

1 *NBTA 


PZE 

NBTA 


BCI 

1 f TR 


PZE 

TR 


BCI 

1 t RHOR 


PZE 

RHOR 


BCI 

1 fFP 


PZE 

FP 


BCI 

1 *FT 


PZE 

FT 


BCI 

1 f ES 


PZE 

ES 


BCI 

1 »TSTAG 


PZE 

TSTAG 


BCI 

1 «RHOS 


PZE 

RHOS 


HTR 

** 

// 

CONTRL 

// 

OMEG 

COMMON 

4500 

F 

COMMON 

30 

CAPM 

COMMON 

30 

A 

COMMON 

300 

G 

COMMON 

900 

SMALE 

COMMON 

900 

CONH 

COMMON 

1 

CONK 

COMMON 

1 

CONR 

COMMON 

1 

CONPRF 

COMMON 

1 

CONNO 

COMMON 

1 

N 

COMMON 

1 

M 

COMMON 

1 

L 

COMMON 

30 

NIT 

COMMON 

1 

EPSl 

COMMON 

1 

FP52 

COMMON 

1 

CAPLAM 

COMMON 

900 

DHFO 

COMMON 

30 

AMC 

COMMON 

1 

ICl 

COMMON 

1 

YSTO 

COMMON 

30 

PI 0 

COMMON 

1 

TlO 

COMMON 

1 

US 

COMMON 

1 

RH02 

COMMON 

1 

EPS5 

COMMON 

1 

NB 

COMMON 

1 

CP 

COMMON 

5 

BETA 

COMMON 

5 
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AM 

COMMON 

5 

TP 

COMMON 

1 

DELT 

COMMON 

1 

IT 

COMMON 

1 

NBTA 

COMMON 

5 

TR 

COMMON 

1 

RHOR 

COMMON 

1 

FP 

COMMON 

1 

FT 

COMMON 

1 

ES 

COMMON 

1 

TSTAG 

COMMON 

1 

RHOS 

COMMON 

END 

1 
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PROGRAM INPUT 


The input necessary to utilize the included program is presented in this 
appendix. The species in the undissociated free stream must he assigned a 
value of the subscript p beginning with unity. Similarly^ the species to be 
considered (atoms, molecules, ions, or electrons) in the dissociated mixture 
must be assigned a value of the subscript i, whereas each elemental particle 
(atomic particle or electron) is given a value of the subscript j. Electrons, 
if considered as a species, must be assigned the value of 1 for the sub- 
script i. The program as it presently stands is limited to the consideration 
of ten elements with five species in the \indissociated free stream and thirty 
species in the dissociated mixture. The capacity of the program with regard 
to the number of species considered is easily increased by changing pertinent 
dimension statements. 

A description of additional input required by the program is given in the 
following table: 
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Input 

Program symbol 

Description 

Unit 

Thermodynamic input 


0MEG 

cth characteristic vibrational frequency of Ith electronic level of 
ith species (<^2c = ^ 

cm"^ 

y± 

YST0 

Initial guess, not equal to zero, for mole number of ith species 




satisfying bj = ^^ij^i 'w^ere bj = product 




P 

apjpp being evaluated in the free stream 


fi 

F 

Zero if species is atom or atomic ion; unity if species is molecule 
or molecular ion 


Mi 

CAPM 

Molecular weight of ith species 

grams/mole 


A 

Number of jth atoms in ith species 


Sii 

G 

Degeneracy of Xth electronic energy level of ith species 


^iX 

SMALE 

Excitation energy of Ith electronic level in ith species 

cm"^ 

h 

C0KH 

Planck' s constant 

erg-sec 

k 

C0KK 

Boltzmann constant 

ergs/^ 

e 

C0NR 

Universal gas constant 

ergs/mole -°K 

Pref 

C0NPRF 

Reference pressure 

dynes/cm^ 

% 

C0M0 

Avogadro's number 

parti cle s / mole 

n 

N 

Number of species i 


m 

M 

Number of types of atoms j appearing in mixture (exclude ions) 


^i 

L 

Number of electronic energy levels considered in ith species 


Nit 

NIT 

Maximum number of iterations in equilibrium subroutine; 
normally 100 



EPSl 

Convergence criterion for absolute test of y^^ in equilibrium 
subroutine; normally 10“T 


^2 

EPS2 

Convergence criterion for relative test on y^^ in equilibrium sub- 
routine; normally €2 = 10 and is not used 


AiX 

CAPIAM 

Product of symmetry number and characteristic rotational tempera- 
ture of Xth electronic level of ith species 
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■ 

MFO 

Standard heat of formation at 0® K of ith species 

ergs/mole 

Ml 

AMC 

Molecular weight of free-stream gas at 300° K 

grams/mole 

c 

ICl 

Maximum number of c's for Zth level of ith species (refers to 
for triatomic species) 


pp 

BETA 

Mole fraction of pth species in free-stream gas at 300° K 


Mp 

AM 

Molecular weight of pth species in free stream 

grams/mole 

£S! 

CELT 

Temperature increment for incident and reflected shock solutions j 
normally 10° K 

°K 

It 

IT 

Option code for additional iteration on temperature in incident and 
reflected shock subroutines j normally 3 


”Pp 

NBTA 

Value of i assigned to pth species in dissociated mixture 


fp 

FP 

Fraction required for incrementing pressure in stagnation point i 

solution; normally 5 X lO"^ 


fp 

FT 

Fraction required for incrementing temperature in stagnation point 
solution; normally 3 X 10"5 


^st 

ES 

Convergence criterion for entropy and enthalpy in stagnation point 
solution; normally 10“^ 



Flow input 


Pi 

PIO 

1 

Free-stream pressure 

atmospheres 

Ti 

TIO 

Free-stream temperature 

°K 

Us 

US 

Incident shock speed 

feet/sec 

P2e 

EH02 

Estimated density behind incident shock 

Nondimenslonal 


EPS5 

Convergence criterion for p in incident and reflected shock 
iterations; normally 10"3 


Tse 

TP 

Estimated temperature behind incident shock 

°K 

Tje 

TR 1 

Estimated temperature behind reflected shock 

°K 

P5e 

RHjiH ^ 

Estimated density behind reflected shock 

Nondimenslonal 


TSTAG 

Estimated temperature behind standing shock 

°K 

P3e. 

BH0S 

1 

Estimated density behind standing shock 

Nondimenslonal 
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Once the input has been detennined for a given initial mixtuxe and species 
to he considered^ duplication for solutions with varying flow conditions, such 
as initial pressure, ten 5 )erature , and shock speed, is not necessary. Computa- 
tions involving these variations require only the respecification of quantities 
given in the foregoing table under Flow input. Similarly, if it is desired to 
determine the effect on the solution of one of the species constants, such as 
the heat of formation of that species, it is necessary to change only that 
value. However, if changes in the mole fractions of the initial mixture are 
made, it is necessary to reestimate values of to insure that the charge 

and mass balance constraints are satisfied. 


Input Loading and Comments 

The input is loaded by a symbolic loader routine, LOAD. Subroutine SYMBOL 
describes the input, giving its symbolic name, nximber of locations required, 
and location in COMMON. Any suitable routine may be used to enter the input by 
appropriately modifying statement 1 in the main program. 

For purposes of output, n + 1 cards are read into the program in front 
of the data to identify the n species in the shock processed mixture. Card 1 
contains the number n (cols. 1 to 3) and cards 2 through n + 1 contain 
alphabetic identification of the species (cols. 1 to 6) - for example, N20, 02, 
and A*^. 

After all the data are loaded, the program reads one card (cols. 1 to 3) 
containing a numeric code for computing option desired: 

0 - incident shock only 

1 - incident and reflected shocks 

2 - incident shock and in-flight stagnation conditions 

3 - incident and reflected shocks, and in-flight stagnation 

conditions 

4 - incident shock, and shock-tube stagnation conditions 

5 - incident and reflected shocks, and shock-tube stagnation 

conditions 

6 - incident shock, shock-tube and in-flight stagnation 

conditions 

7 - incident and reflected shocks, shock-tube and in-flight 

stagnation conditions 

The program also uses a routine MATINV to solve a matrix equation, AX = B, 
where A is a square coefficient matrix and B is a matrix of constant vectors. 
Reference to this routine is found in subroutine ECOM following statement 62. 
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The calling sequence of this routine is shown and hriefly described as follows 
in order to allow replacement by a similar routine, if necessary: 

CALL MATIMV (ABLOCK Ml, BSUM DETERM, IPIVOT, INDEX, 11, O) 

ABLOCK - the first location of two-dimensional array of 
matrix A 

Ml - the location of order of A, 1 ^ Ml ^ 11 

BSUM - the first location of two-dimensional array of con- 
stant vectors B 

I - the n\amber of col-umn vectors 

DETERM - gives value of determinant (not used) 

IPIVOT - temporary storage 
INDEX - temporary storage 

II - the maximum order of A 

0 - a factor used in computing determinant 
At the return to the calling program, x is stored at BSUM. 


Sample Output 

A typical program output for air is given herein, p is in standard atmos- 
pheres, Ug, U]-, U2 are in ft/sec, T is in ^K, and p is nondimensional- 

ized by ^ = 1 . 936 ( 10"“5)£ili§£ = 10"3 Other quantities are nondimension- 

^ cc 

stagnation enthalpies which are in cm^/sec^. 
as follows : 


ft^ 

alized with the exception of the 
The sample computer print-out is 


kl 
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INPLT FCR INCICENT SHCCK 

RHC-1 = C. 1272^3^^8 01 P-1 = C.39SS9999E Cl US= 0.3CC0C0CCE C5 T-1 = b.SCCOCCCOE C3 

BETA 

N2 C02 A 

0.7600CCCCE OC C. 1 ICCCOCCE-OC C. 13CCCCCCE-CC 

ASSLPEC RhC= C.09999999E C2 


CUTPLT 


U2 * 0.276571656 C5 


NC. CF f^AJOR ITbRATICNS = 5 


P 


RFC 


1/Z 


H/ZKT 


S/R 


Ml 


0.9690A850E C3 


FINAL X FBC^ 

C.72S36C7^i6-C2 
C.6273C5C8E CC 
C.369513C2t-C2 
C.11379162E-CC 
C.38751387E-C3 
C.51C65551E-01 
C. 17867831E-C2 
C-727CCei9E-Cl 
C-'^i56<i5799E-C3 
C. 1C506238E-CC 
0.5^371 160E-C3 
C- 1C673^^CE-C3 
C.5563^31CE-C2 
C.AUC5^77E-C3 
0.12752569E-C^ 
C.3A333257E-C2 
C.559CC569fc-C2 
C. 1C381966E-C5 


C. 1629^2696 02 


SPECIES 

E- 

N 

N + 

D 

C + 

C 

C + 

A 

A + 

N2 
N2 + 

02 
NC 
NC + 

C2 

CO 

CN 

CC2 


ITERATICN 


C.5627A6^ie CO C. 671123^66 01 C.A3ll67Clh 02 


C.127761A7E 05 C.3132&C90E 02 


SHCCK TUBE STAGNATION PCINT CATA 

IS CCPPUTEC FROM 

THE FOLLOWING 

STANDING SHCCK DATA 



P 

RFC 

1/Z 

F/ZRT 

S/R 

T 

Ml 

C.1C2237C8E C5 

C.65691168E C2 

C.A867A1AAE-CC 

C.617823C9E 

Cl C.'46117C31b C2 

C.22U9058E 05 

C.31326C9CE C2 


SHCCK TtBE STAGNATICN PCINT ENTHALPY = C.759733A7E 12 


CLTPLT 

PRESStRE = C.11359596E C5 TEPPERATLRE = C.22522939E C5 NO. OF HAJGR ITEkATICjNS = A 

1/Z = C.A86257CAE-CC H/ZRT = C.61796A6AE ?1 RHC = C.93622503E C2 S/R = O.Abll72l3E 02 H = C.75973626H 12 


FINAL X FRCP ITERATION SPECIES 


55171687E-C1 

E- 

6693C937E CC 

N 

38569783E-C1 

N + 

999A1716E-01 

C 

38991915E-C2 

□ + 

A730A658E-C1 

c 

A622861CE-C2 

c+ 

56771533E-C1 

A 

6AA19228E-C2 

A* 

12A83915E-C1 

N2 

1192C3AQE-C2 

N2 + 

1CC52688E-C3 

02 

21869799E-C2 

NC 

A462272At-C3 

NC + 

3A161A65E-C5 

C2 

3023AC67E-C3 

CO 

1251851AE-C2 

CN 

62CC5A85E-C7 

CC2 
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REFLECTEC SHOCK 


ASSLHEC RFC = 0.5SSSSSS9E C2 


CUTPLT 


U-R* O.51850532E 


NC. OF MAJOR 

ITERATIONS = 

5 




P 

RFC 

l/Z 

H/ZRT 


S/R 

T 

Ml 

0.1A53986CE 05 

C-IC319477E C3 

C.^712C168E-CC 

C.62550C5^6 

Cl 

0.A77CC31AE C2 

0.253AA&2TE C5 

V..31326C9CE C? 


FINAL X FROM ITERATION 

SPECIES 

C.77387C73E-C1 

£- 

0.6ii283A80t CC 

N 

0.556C6787E-C1 

N + 

C.95394575E-C1 

0 

0.613075CSE-C2 

0 + 

C.-^5581391E-C1 

C 

0.5^3657C6E-C2 

c+ 

C.52A83761E-C1 

A 

0.87725917E-C2 

A + 

C.t63^82l7E-C2 

N2 

C.IC2772A7E-C2 

N2 + 

C.7822S235E-CA 

02 

C.1AC52661E-C2 

NO 

0.^1299297b-C3 

NO* 

C. 17578935E-C5 

C2 

0.16513C67E-C3 

CO 

0.6A58AA81E-C3 

CN 

C.23657988E-C7 

CC2 


IN FLIGHT 
STAGNAT ICN 


STAGNATION 
ENTHALPY = 


POINT DATA IS CCHPUTEC FRC^ THE INCinENT 
C.‘'.Cfcg66eit 12 stagnation LNTRLPY = C 


SHOCK CLTPUT 
^T1167Cll- C2 


CtTPLT 

PRESSORE = C.lCligiAAE TEMPERATURE = C.12a5A92gt C5 

l/Z = C.5617g306t CC H/ZRT = 0.67CC527AE Cl RHO = C.16882119E Q2 


NC. OF MAJOR ITERATILNS = 


S/R = ?.A311A19^E 02 


U.^C69S12'«£ 12 


L X FRCM ITERATION 

SPEC 

7^532A3AE-C2 

E- 

629CC783E CC 

N 

3807A829E-C2 

N + 

1136^757E-CC 

0 

39727895E-G3 

0 + 

51GAC16CE-C1 

C 

18QA^C82£-C2 

C + 

7256C759E-01 

A 

^7259831E-C3 

A + 

,ICA22159E-CC 

N2 

55593C66E-C3 

N2 + 

.10777555E-C3 

U2 

.55677C62C-C2 

NO 

A1558252E-C3 

NQ + 

12763625E-C^ 

C2 

,33A9?790t-02 

CO 

5576579AE-G2 

CN 

1C233709E-05 

C02 
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APPENDIX C 


PHYSICAL CONSTANTS 


The physical constants required hy the computer program for use with the 


27 chemical species are as follows: 

Planck* s constant^ h^ erg-sec 6.62517 X 10’^^ 

Boltzmann constant^ ergs/^K 1.38o44 X 10“1^ 

Universal gas constant^ fly ergs/mole-^K 8.31469 X 10^ 

Reference pressure^ Pref^ dynes/cm^ I.OI325 X 10^ 

Avogadro*s number^ particles/mole 6.02322 X 10^5 

Speed of lights c^ cm/sec 2.99793 X lO^O 


The molecular weight and heat of formation of the chemical species are 
given in the following table: 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


Ml 


0 

4.70729 

18.72607 

47.28829 

2.46741 

15.60589 

49.47996 

1.05410 

7.11258 

17.97182 

41.49225 

5.89944 

0 

15.20255 

41.85170 

0 

15.05556 

0 

11.62808 

-.96252 

.89860 

9.82405 

-1.15815 

12.58567 

4.56056 

- 5.95146 

.84975 


10I2 


If electrons are included as a species^ they must be assigned a value of 
1 for the subscript i; otherwise values of i may be assigned to species as 
desired. 
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APPENDIX C 


The number of jth particles in ith species is shown in the following table 


Species 

i 



for - 


N 

0 

c 

A 

e” 

e" 

1* 

0 

0 

0 

0 

+1 

N 


1 

0 

0 

0 

0 

N+ 


1 

0 

0 

0 

-1 

N++ 


1 

0 

0 

0 

-2 

0 


0 

1 

0 

0 

0 

0+ 


0 

1 

0 

0 

-1 

0++ 


0 

1 

0 

0 

-2 

0- 


0 

1 

0 

0 

+1 

c 


0 

0 

1 

0 

0 

c+ 


0 

0 

1 

0 

-1 

C++ 


0 

0 

1 

0 

-2 

c- 


0 

0 

1 

0 

+1 

A 


0 

0 

0 

1 

0 

A+ 


0 

0 

0 

1 

-1 

A++ 


0 

0 

0 

1 

~2 

Ng 


2 

0 

0 

0 

0 



2 

0 

0 

0 

-1 

Og 


0 

2 

0 

0 

0 

0g+ 


0 

2 

0 

0 

-1 

Og" 


0 

2 

0 

0 

+1 

NO 


1 

1 1 

0 

0 

0 

N0+ 


1 

1 

0 

0 

-1 

CO 


0 1 

1 

1 

0 

0 

C0+ 


0 1 

1 

1 

0 

-1 

CN 


1 

0 

1 

0 

0 

COg 


0 

2 

1 

0 

0 

NgO 


2 

1 

0 

0 

0 


1’Values of j for the elements shown may be assigned as desired. 

If electrons are included as a species^ they must be assigned a value of 
1 for the subscript i; otherwise values of i may be assigned to species as 
desired. 
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Spectroscopic constants and quantities for the 27 species are as follovs 
(values of ^Iq, ^±l "being listed only where required): 


Species 

i 

1 

Siz 

Hi 


Species 

i 

1 

Hi 

HI 


Species 

i 

1 

Hi 

SiZ 

e- 

1* 

1 

2 

0 


N+ 


1 

9 

0 


N++ 


1 

6 

0 









2 

5 

15316 




2 

12 

57280 

N 


1 

4 

0 




3 

1 

52687 




3 

10 

101026 



2 

10 

19225 




4 

5 

47168 




4 

2 

151044 



5 

6 

28840 




5 

15 

92245 




5 

6 

145920 



if 

12 

85550 




6 

9 

109220 




6 

4 

186802 



5 

6 

86180 




7 

5 

144189 




7 

10 

205078 



6 

12 

88 i 4 o 




8 

12 

149000 




8 

2 

221502 



7 

2 

93582 




9 

3 

155130 




9 

6 

250507 



8 

20 

94800 




10 

3 

164612 




10 

6 

245690 



9 

12 

95500 




11 

18 

166650 




11 

10 

267242 



10 

4 

96752 




12 

3 

168895 




12 

12 

287650 



11 

10 

96810 




13 

9 

170620 




13 

6 

297210 



12 

6 

97800 




14 

5 

174212 




14 

2 

501088 



13 

10 

99660 




15 

1 

178274 




15 

32 

510500 



lif 

18 

104000 




16 

4 i 

187090 




16 

4 

514224 



15 

34 

104700 




17 

16 

189100 




17 

22 

317550 



;i6 

56 

105000 




18 

3 

190121 




18 

14 

320288 



^17 

22 

106600 




19 

12 

197200 




19 

10 

521040 



l8 

16 

107200 


! 


20 

30 

205000 




20 

30 

529000 



19 

l6 

107600 




21 

21 

205700 




21 

22 

333500 



20 

54 

109600 




22 

153 

210500 




22 

22 

356290 



21 

12 

109860 




23 

12 

214828 




23 

20 

339550 



22 

106 

110550 




24 

12 

218000 




24 

48 

345000 



23 

132 

111600 




25 

546 

220500 




25 

24 

548000 



24 

136 

112200 




26 

37 

225000 




26 

32 

564000 



25 

6 

112510 




27 

15 

226000 




27 

12 

568620 



26 

68 

112820 




28 

15 

228000 




28 

22 

574000 



27 

30 

112910 




29 

i 4 

250500 




29 

52 

579000 



28 

2 

115600 




30 

5 

254250 




30 

2 

580700 



29 

10 

115900 















30 

l 4 

117000 














•¥r 

If electrons are included as a species^ they must he assigned a value of 
1 for the subscript i; otherwise values of i may be assigned to species as 
desired. 
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Species 

i 

Z 

^iZ 


0 


1 

5 

0 



2 

3 

159 



3 

1 

227 



4 

5 

15868 



5 

1 

33792 



6 

5 

73768 



7 

3 

76795 



8 

15 

86630 



9 

9 

88630 



10 

5 

95476 



11 

3 

96226 



12 

36 

97420 



13 

15 

97488 



l 4 

15 

99094 



15 

9 

99680 



16 

9 

100000 



17 

15 

101140 



18 

5 

102116 



19 

3 

102412 



20 

5 

102662 



21 

25 

102865 



22 

15 

102908 



23 

9 

103869 



24 

15 

104000 



25 

25 

105385 



26 

15 

105408 



27 

56 

106000 



28 

56 

107000 


Species 

i 

Z 

Sil 

Hi 

O'*' 


1 

4 

0 



2 

10 

26820 



3 

6 

4 o 465 



4 

12 

120000 



5 

10 

165990 



6 

12 

185400 



7 

6 

189000 



8 

2 

195710 



9 

2 

203942 



10 

42 

207600 



11 

26 

212800 



12 

2 

226851 



13 

24 

230000 



l 4 

80 

232700 



15 

16 

233900 



16 

18 

239600 



17 

4 o 

245500 



18 

20 

248000 



19 

6 

250300 


1 

20 

32 

251900 



21 

44 

254000 



22 

84 

255500 



23 

116 

256000 



24 

18 

258100 



25 

10 

259300 

1 


26 

72 

261500 



27 

166 

265500 



28 

134 

270000 



29 

142 

276500 



30 

24 

283000 


Species 

i 

Z 

SiZ 


0++ 


1 

9 

0 



2 

5 

20271 



3 

1 

43184 



4 

5 

60312 



5 

15 

120050 



6 

9 

142383 



7 

5 

187049 



8 

3 

197087 



9 

3 

210459 



10 

12 

270000 



11 

9 

283900 



12 

26 

294000 



13 

14 

303000 



l 4 

1 

313801 



15 

50 

327000 



16 

10 

332000 




15 

338700 



18 

1 

343303 



19 

9 

350200 



20 

12 

357500 



21 

46 

565000 



22 

48 

570500 



23 

158 

380000 



24 

62 

394000 



25 

90 

398000 



26 

558 

403400 



27 

72 

425000 



28 

311 

430000 



29 

48 

438000 

, 


30 

15 

442710 
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Species 

i 

1 

Siz 

Hi 

0“ 


1 

6 

0 

c 


1 

9 

0 



2 

5 

10194 



5 

1 

21648 



4 

5 

33735 



5 

9 

60360 



6 

3 

61982 



7 

15 

64090 



8 

3 

68858 



9 

15 

69700 



10 

3 

70744 



11 

9 

71365 



12 

5 

72611 



13 

1 

73976 



l 4 

9 

75256 



15 

5 

77681 



l6 

9 

78130 



17 

21 

78230 



18 

18 

78320 



19 

10 

78600 



20 

9 

79318 



21 

18 

8 o 4 oo 



22 

12 

81200 



23 

5 

81770 



24 

1 

82252 



25 

53 

83800 



26 

103 

84000 



27 

30 

84940 



28 

6 

85400 



29 

503 

86400 



30 

19 

86500 


Species 

i 

1 

SlZ 


c+ 


1 

6 

0 



2 

12 

43030 



3 

10 

74931 



4 

2 

96494 



5 

6 

101800 



6 

6 

110650 



7 

2 

114900 



8 

2 

116558 



9 

l 4 

119400 



10 

6 

151751 



11 

12 

135800 



12 

20 

136000 



13 

4 

142024 



l 4 

10 

145551 



15 

10 

150465 



16 

2 

157234 



17 

6 

162522 



18 

12 

167000 



19 

10 

168124 



20 

20 

168900 



21 

2 

173548 



22 

6 

175293 



23 

16 

178350 


1 

24 

32 

179000 



25 

26 

182000 



26 

4 

184689 



27 

12 

186450 



28 

10 

188600 



29 

30 

195500 



30 

20 

196572 


Species 

i 

1 

giz 

HI 

C++ 


1 

1 

0 



2 

9 

52360 



5 

3 

102351 



4 

9 

137420 



5 

5 

145875 



6 

1 

182520 



7 

3 

238161 



8 

1 

247170 



9 

3 

258951 



10 

9 

259662 



11 

15 

269960 



12 

5 

276843 



13 

15 

309100 



l 4 

1 

311721 



15 

12 

318700 



16 

61 

522550 



17 

5 

324212 



18 

12 

528000 



19 

31 

535000 



20 

15 

337626 



21 

13 

340000 



22 

7 

341368 



23 

3 

343256 



24 

30 

545000 



25 

60 

346600 



26 

7 

348000 



27 

12 

376600 



28 

18 

381500 



29 

13 

384545 



30 

34 

386000 
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Species 

i 

1 

SiZ 

Hi 

A+ 


1 

6 

0 



2 

2 

108723 



5 

20 

152400 



4 

12 

134800 



5 

6 

138600 



6 

28 

142700 



7 

6 

145200 



8 

12 

147650 



9 

10 

148750 



10 

10 

150000 



11 

12 

155160 



12 

50 

158500 



15 

12 

161000 



l 4 

2 

167309 



15 

l 4 

170600 



l 6 

26 

175000 



1 ? 

38 

174800 



18 

6 

179700 



19 

4 o 

185000 



20 

6 o 

186000 



21 

62 

190200 



22 

54 

192200 



23 

120 

194000 



24 

84 

196000 



25 

48 

200000 



26 

18 

205000 



27 

6 

208595 



28 

220 

210000 



29 

50 

215000 



50 

70 

220000 


Species 

i 

z 

SiZ 

Hi 

A++ 


1 

5 

0 



2 

5 

inp 



5 

1 

1570 



4 

5 

i 4 oio 



5 

1 

35267 



6 

9 

114400 



7 

10 

128000 



8 

28 

144650 



9 

15 

156950 



10 

5 

174375 



11 

15 

182000 



12 

75 

189500 



13 

15 

196600 



14 

5 

200000 



15 

18 

204700 



16 

21 

208500 



17 

27 

210800 



18 

53 

214500 



19 

46 

224500 



20 

9 

231500 



21 

5 ! 

235000 



22 

26 

240000 



25 

50 

246056 



24 

91 

251000 



25 

24 

258000 



26 

127 

270000 



27 

121 

278500 



28 

116 

282000 



29 

111 

286000 



30 

524 

302000 


Species 

i 

z 

siz 

Hi 

c- 


1 

6 

0 

A 


1 

1 

0 



2 

5 

93144 



3 

3 

93751 



4 

1 

94554 



5 

3 

95400 



6 

5 

104102 



7 

12 

105500 



8 

8 

106150 



9 

1 

107054 



10 

8 

107220 



11 

4 

108000 



12 

4 

111280 



13 

16 

111750 



14 

8 

112900 



15 

20 

113550 


1 

16 

16 

114750 



17 

8 

115000 



18 

3 

116660 


1 

19 

12 

116960 



20 

8 

117170 



21 

1 

117563 



22 

16 

118530 



23 

36 

119300 



24 

56 

120250 


i 

25 

20 

120700 



26 

52 

121750 



27 

36 

122200 



28 

56 

122700 



29 

20 

125500 



30 

28 

124136 
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Species 

i 

1 

Siz 

^iZ 


^Z 1 

N 2 


1 

1 

0 

5.725 

2343.9 



2 

3 

49757 

4.125 

1446.5 



5 

6 

5931 *^ 

4.687 

1719.6 



k 

2 

68953 

4.630 

1625.9 



5 

1 

70700 

4.235 

1518.0 



6 

6 

87984 

5.226 

2018.0 



1 

2 

0 

5.531 

2191.0 



2 

4 

9020 

4.929 

1887.9 



3 

2 

25570 

5.966 

2396.7 



4 

2 

64550 

4.676 

2035.1 

O2 


1 

3 

0 

4.137 

1568.3 



2 

2 

7882 

4.080 

1496.4 



3 

1 

15121 

4.004 

i 4 i 8.7 



4 

3 

35713 

2.970 

796.5 



5 

1 

36213 

2.347 

633.3 



6 

3 

49363 

2.341 

692.4 

O2'" 


1 

4 

0 

4.783 

1859.9 



2 

8 

31500 

3.156 

1025.3 



3 

4 

38300 

3.028 

886.6 



4 

4 

48100 

3.673 

1179.7 

02’ 


1 

4 

0 

3.430 

1286.0 



2 

4 

13400 

2.767 

975.0 



3 

4 

24200 

2.609 

547.0 

NO 


1 

2 

0 

2.440 

1890.1 



2 

2 

121 

2.440 

1889.7 



3 

2 

44200 

2.860 

2358.3 



4 

4 

45440 

1.609 

1030.1 



5 

2 

53290 

2.866 

2301.0 



6 

4 

52376 

2.859 

2380.0 



7 

2 

60860 

2.845 

2357.8 



8 

4 

60020 

1.900 

1200.7 

NO"^ 


1 

1 

0 

2.866 

2360.8 



2 

6 

39982 

2 . 4 o 4 

1725.5 



3 

3 

58523 

1.902 

1210.5 



4 

6 

72384 

1.801 

1132.4 



5 

2 

73084 

2.266 

1585.6 

CO 


1 

1 

0 

2.766 

2156.8 



2 

6 

48474 

2.405 

1724.8 



3 

3 

55380 

1.904 

1208.5 



4 

6 

61785 

1.803 

1130.2 



5 

2 

64747 

2.303 

1498.4 



6 

3 

83831 

2.962 

2184.5 



7 

1 

86918 

2.802 

2070.0 

CO+ 


1 

2 

0 

2.831 

2199.1 



2 

4 

2 o 4 o 7 

2.273 

1548.5 



3 

2 

45634 

2.568 

1706.3 

CN 


1 

2 

0 

2.7207 

2055.56 



2 

4 

9115 

2.4571 

1801.55 



3 

2 

25798 

2.8186 

2143.88 

CO2 


1 

1 

0 

1.121 

667.3 
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(a) Traveling normal shock. 
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(h) Reflected normal shock. 




Figure l.~ Flow configurations. 
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(a) Pressure. 

Figure 2.- Thermodynamic ratios across traveling normal shock for carton and argon 

free air. T. = 500° K. 
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(a) Pressure. 

Figure Shock- tube stagnation -point conditions for carbon and argon free air 
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(b) Temperature. 


Figure 3.- Continued. 
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Figure 4.- Electron number density behind incident shock in particles per cubic centimeter 

for carbon and argon free air. T-i = 300° K. 



Figure 5-“ Molecular wei^t 'behind incident shock for carhon and argon free air. 

Ti = 300° K. 
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(c) Density. 
Figure 6.- Concluded. 
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Pressure ratio, pg/p^ 



(a) Press\ire. 

Figure 7 *~ Thermodynamic ratios across refiected. normal shock for carbon and. argon free air. 

Ti = 273-2° K. 
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Incident shock speed, U^, m/s 

(c) Density. 

Figure 8.- Concluded. 
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"The aeronautical and space activities of the United States shall he 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof." 
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